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ABSTRACT 

We develop and analyze a simple cellular automaton (CA) model that re- 
produces the main properties of the evolution of soft X-ray coronal loops. We 
are motivated by the observation that these loops evolve in three distinguishable 
phases that suggest the development, maintainance, and decay of a self-organized 
system. The model is based on the idea that loops are made of elemental strands 
that are heated by the relaxation of magnetic stress in the form of nanofiares. 
In this vision, usually called "the Parker conjecture" (Parker 1988), the origin 
of stress is the displacement of the strand footpoints due to photospheric con- 
vective motions. Modeling the response and evolution of the plasma we obtain 
synthetic light curves that have the same characteristic properties (intensity, 
fluctuations, and timescales) as the observed cases. We study the dependence of 
these properties on the model parameters and find scaling laws that can be used 
as observational predictions of the model. We discuss the implications of our 
results for the interpretation of recent loop observations in different wavelengths. 

Subject headings: Sun: corona - Sun: flares - Sun: magnetic topology - Sun: 
X-rays, gamma rays 

1. Introduction 

EUV and soft X-ray observations of the solar corona reveal a high degree of structuring 
and intermittency, especially in and around active regions (ARs). Due to the frozen-in 
condition of the coronal plasma, the emitting material is ordered by the magnetic field 
in elongated features commonly known as coronal loops. The study of the structure and 
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evolution of loops is therefore fundamental for understanding coronal dynamics, particularly 
in relation to the origin of coronal heating (Mandrini et al. 2000). Historically, the first 
attempts to physically describe coronal loops began with what is known as the Rosner, Tucker 
& Vaianna (RTV, 1978) approach, which considers loops as monolithic structures containing 
plasma in static equilibrium. Although soft X-ray loop observations seem to be consistent 
with the RTV assumptions (Porter & Klimchuk 1995), more recent observations from the 
Transition Region and Coronal Explorer (TRACE) suggest that EUV loops are actually too 
dense to be in static or quasi-static equihbrium (Aschwanden et al. 2001, Winebarger et 
al. 2003) . The properties of these loops are more consistent with the predictions of models 
based on impulsive heating (Winebarger et al. 2003, Klimchuk 2006). 

TRACE observations also show a high degree of filamentation of the coronal structure, 
suggesting that what appear to be monolithic loops can actually be collections of unre- 
solved individual strands with more or less independent evolutions. Although the question 
of whether coronal loops are isothermal or multithermal has been a matter of a heated 
debate during several years (Aschwanden & Nightingale 2005, Schmelz & Martens 2006), 
recent observations seem to indicate that both kinds of loops actually occur (Schmelz et al. 
2009). Combining all the available information, there are strong arguments in favor of many 
loops being made of unresolved strands that evolve more dynamically than what quasi-static 
models predict (for a recent review on the subject see Khmchuk 2009). 

There are two main families of mechanisms proposed to explain the origin of coronal 
heating: the so called AC and DC heatings. Generally speaking, AC refers to the dissipation 
of MHD waves, while DC relates to the dissipation of highly stressed magnetic fields. The 
latter involve the release of free magnetic energy at locations of strong gradients in the 
magnetic field, i.e., current sheets. Ultimately, the source of energy is the distortion of the 
magnetic structures by convective motions. This is precisely the mechanism envisioned by 
Parker (1988). He proposed that, since the footpoints of magnetic strands are being displaced 
by photospheric granular motions, there must be a continuous increase of the magnetic stress 
between neighbor strands. These stresses accumulate until some threshold is reached and 
reconncction occurs. In the process, stored magnetic energy is released and the plasma is 
heated. Because of similarities with the usual flaring process and the energy scale involved, 
these reconnection events are usually called nanoflares (Cargill 1994, Cargill & Khmchuk 
2004). One key point of the process described above is the existence of a threshold or critical 
value for the reconnection to occur. A theoretical support for the existence of this threshold 
is given in Dahlburg et al. (2005, 2009). They found that when certain field misahgnment is 
reached the secondary instability provides the mechanism for a rapid energy release. 

The intermittent and filamentary behavior, and the possible presence of power law 
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distributions (see e.g., Aschwanden & Parnell 2002), led some authors to propose that the 
solar corona might be a good example of a self-organized critical (SOC) system (Bak 1990). 
Beginning with the seminal work of Lu & Hamilton (1991), there have been an important 
number of studies involving the use of SOC models in relation to the coronal heating problem 
(for a review on SOC models see Charbonneau et al. 2001). 

In our recent paper, Lopez Fuentes et al. (2007), we studied the evolution of coronal 
loops observed with the Solar X-ray Imager (SXI) on board the Geosynchronous Opera- 
tional Environmental Satellite 12 (GOES-12). Since this instrument observes the sun almost 
continuously, it allowed us to follow the evolution of a set of coronal loops from birth to disap- 
pearance. We found that the loop evolution can be separated into three main phases, which 
we called: rise, main and decay. As we indicated there and began to explore in Klimchuk, 
Lopez Fuentes & DeVore (2006), the observed evolution can be related to the development, 
maintenance and destruction of a critical system. Here, we develop a simple model based on 
a cellular automaton (CA) that reproduces the observed evolution. We model the plasma 
response to the heating using the EBTEL hydrodynamics code (Klimchuk et al. 2008). We 
then create synthetic light curves taking into account the temperature-dependent response of 
the SXI instrument. We compare these with the observed light curves and analyze how the 
different parameters of the model influence the properties of the obtained light curves, such 
as the characteristic intensities and timescales, which we studied in detail in Lopez Fuentes 
et al. (2007) for the observed cases. 

In Section[2]we sumarize the results from Lopez Fuentes et al. (2007) and show examples 
of observed light curves that will be later compared with the model. In Section |3] we describe 
the CA model and in Section H] we analyze our results. We discuss and conclude in Section [51 

2. Observed loop evolution 

In Lopez Fuentes et al. (2007) we studied the evolution of a set of coronal loops observed 
with the Solar X-ray Imager (SXI) on board GOES-12 (Hill et al. 2005). This instrument 
has the advantage of observing the sun almost continuously, allowing us to follow the whole 
evolution of loops. Previous similar studies based on Yohkoh/SXT or TRACE observations 
were limited by the occultation times of the satellites and times of interrupted observation due 
to high radiation levels (e.g., when passing through the South Atlantic anomaly), therefore 
reducing the continuous observation times to a small portion of the loop lifetime (see e.g. 
Porter & Klimchuk 1995). 

We analyzed a set of 457 images taken with the telescope in the filter position "open" , 
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which is most sensitive to emission from plasmas below 2MK. The analyzed data span three 
days (August 26-28, 2003) of observations with a cadence of approximately 6 images per 
hour. From these data we selected and tracked the evolution of 7 coronal loops. Following 
a procedure that is thoroughly described in Lopez Fuentes et al. (2007), we obtained light 
curves of the intrinsic (background subtracted) intensity of each of these loops. We found 
that the loop evolution can be separated into three parts: a rise phase during which the loop 
intensity slowly increases, a main phase during which the intensity remains approximately 
constant (or its long term variation is much less pronounced), and a decay phase during 
which the intensity decreases until the loop is almost indistinguishable from the background. 
In Figure [1] (upper panels) we reproduce adapted versions of two of the light curves studied 
in Lopez Fuentes et al. (2007). In the original analysis we applied a procedure to remove the 
background contribution that left a residual component of long term variation (see discussion 
in Section 3.1 in Lopez Fuentes et al. 2007). For an easier comparison with the model light 
curves, here we removed the residual background by subtracting in each case a linear function 
of time, so that the initial and final intensities of the loops are reduced nearly to zero. 

We found that the durations and characteristic timescales of the rise and decay phases 
are longer than the expected cooling times. We concluded that this is consistent with two 
alternate scenarios. In the first one, loops are monolithic and the heating source varies very 
slowly with time, so the plasma responds quasi-statically. In this way, the three evolution 
phases would reflect the phases of the long term change of a slowly varying heating source. 
In the second scenario, loops are made of unresolved and independent magnetic strands 
which are heated by impulsive events such as nanoflares. Since the observed loop emission is 
the summed contribution of the constituent strands, the different phases would be directly 
related to the variation of the number and intensity of the energy dissipation events that 
heat the strands. We argued that under this particular scenario the observed loop evolution 
can be related to the development, maintenance and destruction of a self-organized system. 
In the next sections we present and analyze a model that reproduces precisely this kind of 
behavior. 



3. Cellular automaton model 

3.1. General description 

In this Section we describe the construction of a cellular automaton (CA) model that 
reproduces the main features of the observed loop evolution. Motivated by the arguments 
given in Section [H the model simulates the process described there regarding the scheme: 
motion of magnetic strand footpoints — ?■ increase of magnetic stress — )■ energy release by 
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reconnection. Figure 2 shows a cartoon that illustrates the main idea. 

Let us begin with a bunch of thin parallel magnetic flux tubes or strands that connect 
two different regions of the solar surface (panel (a)). For simphcity the strands are drawn 
straight, so the top and bottom planes correspond to two relatively distant photospheric 
regions of opposite magnetic polarity. As time goes on, strand footpoints are horizontally 
dragged by random photospheric motions, so the strands become progressively more mis- 
aligned. This produces growing magnetic stresses and intensifying current sheets in the region 
separating adjacent strands (panel (b)). For two given neighbor strands, the misalignment 
angle increases until a critical value is reached (i.e., the secondary instability develops) and 
the energy is released in the form of an impulsive event. Consequently, the field between the 
two strands relaxes. This sudden change in the field may occasionally cause new instability 
conditions between these strands and their other neighbors. Under some conditions instabil- 
ities propagate farther producing avalanches of events involving many strands. This is the 
typical situation considered in SOC models. For clarity, we adopt the following terminology. 
The basic interaction between two strands that occurs each time the critical value is reached 
is called a "reconnection event." The total release of energy into a single strand by its re- 
connection with one or more neighbors over a short period of time is called a "nanofiare." 
This follows the convention used by the loop modeling community. 

The system continues to evolve in response to the footpoint driving until there is a 
rather abrupt appearance of a critical state in which the energy input by the photospheric 
motions is statistically balanced by the dissipation due to reconnection events. In the critical 
state the total energy of the system fiuctuates intermittently around a constant level. Local 
intermittency results from the occurrence of nanofiares of different sizes. 

The idea is then to develop a numerical algorithm that quantitatively reproduces the 
above scheme. First of all, we need to translate footpoint displacements to magnetic field 
stresses. To accomplish this we begin with the simple scenario shown in Figure 2 (c). Let us 
concentrate on one particular strand and, for simplicity, suppose that only one footpoint (the 
lower one) is being displaced. If the particular strand has a vertical magnetic field strength 
By, any horizontal displacement of its footpoint will produce a new horizontal component of 
the field 

5Bh^ Butane, (1) 

where 6 is the angle that the now inclined strand forms with the vertical direction. It is easy 
to see that if the strand has a length L and the footpoint has been displaced a distance d, 
then tan^ = d/L. Therefore: 
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5B^, = BJ/L. (2) 

The distance d is a characteristic length scale of the photospheric motions, which can be 
associated with the sustained displacement of a strand footpoint during the turnover time 
of a convective granule (~ 10^ s). Assuming a typical photospheric velocity v ^ 1 km/sec, 
results in d ^ 1000 km. After a number of successive displacements the strand would acquire 
a horizontal field 

Bh = T.5Br,. (3) 

At this point we have assumed that B^ increases with each and every step. This is only true 
if the new step is in the same direction as the previous step or if the change in direction 
causes the strand to wrap around a neighbor. The latter will generally be true if the random 
walk step size is larger than the characteristic separation of neighbor footpoints. However, 
even then, some steps will cause the strand to back track and B^ will decrease. We treat 
this more realistic situation below. 

As a consequence of the random walk, the horizontal field B^ will gradually increase 
until the critical condition referred to above is reached. If we call 9c the critical angle of 
misalignment between two adjacent strands, it is easy to see that there is an associated 
critical value of the difference in their horizontal field components: 

Bc = B^.tanOc. (4) 

Note that 6*^ refers to the relative angle between the adjacent strands and not the inclination 
from vertical, 6, in Equation 1. 

Assuming ID displacements (see Section I3l2|) . we relate the misalignment between two 
neighbor strands (with indices p and q) to the difference of their horizontal field components 
ABh = Bh{p) — Bh{q). The fulfillment of the critical condition implies |Ai?/i| > B^.. If we 
assume that all the field is equally distributed after reconnection, the new (primed) horizontal 
field in each strand will respectively be: 

B'^ip) = B^{p) ~ AB,,/2 and B'^{q) = Bf,{q) + ABj,/2. (5) 

The released energy corresponds then to the difference in the total magnetic energy 
before and after reconnection. It can be easily demonstrated that 



7 




1 \ {AB,f 
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per unit volume, where we have assumed that the energy is divided equally between the 
plasma in the two strands. 



For the numerical implementation of the model we use a ID array of points. Each 
point corresponds to a single magnetic strand and, initially, all strands have horizontal 
magnetic field = 0. At each time step, a small random amount of field 6Bh is added to 
each site according to a statistical distribution that simulates the photospheric displacements 
of strand footpoints. Although the field increase in the model is ID, we apply a distribution 
based on 2D footpoint displacements. The idea is illustrated in Figure [3]-a. If one compares 
the position of a strand footpoint at time step (t — 1) with its position at time step t (gray 
arrow), the probability that during the following time step, (black arrows) the footpoint 

moves farther away increasing the horizontal magnetic strength is approximately 3/4, while 
the probability that it moves back on its previous track, canceling (or almost canceling) 
the previous displacement is roughly 1/4. Therefore, we add the SB^ values according to 
the distribution shown in Figure [3]-b. Values from the gaussian distribution centered on 
are 3 times more probable than values given by the gaussian centered on —do. The distance 
do = 1000 km is the average horizontal displacement of a strand footpoint during the turnover 
time of a granular cell (in our scheme, the duration of each time step). The width of the 
gaussians is 20% of do- As explained in the previous Subsection, the random displacements 
taken from this distribution translate into magnetic inputs through Equation [2J 

It is easy to see that under this simple ID scheme all strands will eventually acquire 
positive values of B^. To keep the simplicity of the model and to be consistent with the fact 
that neighbor strand footpoints tend to move in different directions, we alternate the sign of 
the displacements d among neighbor strands. This means that, identifying each mesh site 
with an index i, we add the random values obtained from the distribution to strands with 
even i, and we subtract them from strands with odd i. Thus, in the long run even strands 
acquire positive B^ while odd strands acquire negative B^. This alternation of SB^ simulates 
in ID the magnetic stress produced by footpoints moving in different random directions. 

The efficiency of the process that injects "horizontal component" of the magnetic field 
through the motion of the footpoints depends on the average distance between neighboring 
strands. If the strands are closer than c/q, then the process can be considered efficient in 



3.2. Model implementation 



- 8 - 



tangling the strands and accumulating magnetic stress. However, if there is no injection of 
new flux, then the distance between strands will increase as the footpoints disperse. The 
driving process will become progressively less efficient. To simulate the decay of the driving 
mechanism due to footpoint dispersion, at a chosen point during the system evolution we 
make the distribution shown in Figure |3]-b slowly evolve to the symmetrical distribution 
shown in Figure [3]-c. In the final state, the driver distribution is completely inefficient at 
creating new Bh. As we will see below, this evolution of the driver leads to the decay of 
the system. The start time and the duration of the evolution from the initial to the final 
distribution are set as parameters of the model. 

In Subsection 13.11 we defined Be = B^tanOc. Therefore, tan 6'c is one of the relevant 
parameters of the model. At each time step, after adding the random SB^ values to all the 
strands, we compare Bh on each strand with its immediate neighbors. For a strand having 
index i we compute the difference ABh{i) = Bh{i) — Bh{i + 1). We do this for all the strands 
in the array. We use periodic boundary conditions, so strand i = Nm is compared with 
strand i = I, and in that case ABh{Nm) = Bh{Nm) — -8^(1). Next, we identify all the indices 
where \ABh{i)\ > Be- For all the sites where this happens, we redistribute the horizontal 
field according to Equation O 

As we discussed in Subsection 13.11 each field redistribution is associated with a recon- 
nection event that frees an energy AE, in erg cm~^, according to Equation [61 This energy 
is used to heat the plasma in both strands involved in the event. Before proceeding to the 
next time step we repeat the previous test and corresponding field redistribution as many 
times as necessary until all \ ABh{i)\ values are less than Be. Once this state is reached, the 
system is allowed to evolve to the next time step, a new set of random SBh values is added 
to the Bh of the strands, and the process described above is repeated. 

During the evolution of the system we keep record of the energy dissipated in each 
strand at each time step. Here, we assume that the reconnection events (AE) occur on a 
timescale that is much shorter than the time step duration 10'^ s). As we defined in 
Section 13.11 the sum of all the AE "energy packages" that heat the plasma in the strand 
during a given time step is called a nanofiare. The final output of the CA model is the array 
of nanofiare energies dissipated in each strand as a function of time. 

In the above description we assumed that when the critical condition (|Ai?ft| > Be) 
is fulfilled all the available field ABh is distributed (reconnected) among the two strands. 
This is not realistic. On the Sun, each strand of the tangled field is wrapped around many 
other strands. It is only the horizontal field component along the section of mutual contact 
that is free to reconnect and heat the plasma. We account for this in our simple model 
by introducing a parameter /. The initial response to the reconnection is to change the 
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horizontal field locally by an amount ABh/'i. This change is subsequently spread out along 
the entire length of strand. The adjustment happens at the Alfven speed and is very rapid. 
Assuming a minimum Alfven speed of 1000 km for the solar corona gives a time of 100 s 
or less. The end result is that the horizontal field changes by fABh/2 along the entire strand. 
It can be demonstrated that in this case the energy release on each of the two strands per 
unit volume is 



once again reach the critical condition through continued driving. If / is small, the strands 
will never be far from critical and nanoflares will recur frequently. If / is large, it will take 
many steps to once again reach the critical level of stress. Of course the nanoflares will be 
weaker in the former case than in the latter, so the temporally averaged heating will not 
be greatly different. There are small differences though, reflected in the / dependence of 
equation (fTT!) discussed later. 

For the analysis, we vary the values of the model parameters B^, L, tan6'c, /, N^, and 
r. The parameter r is the duration of the nanoflares as defined later in Section 13.31 Table [1] 
shows the different values used. For the vertical magnetic field we use values in the range 
20-120 Gauss, typical of AR loops (Mandrini et al. 2000). We vary the length of the loops 
also considering typical AR values (20-120 Mm). The tangent of the critical angle (tan^^^) 
varies between 0.1 and 1.0, corresponding approximately to angles: 6 deg, 11 deg, 22 deg, 31 
deg, 39 deg, and 45 deg. The reconnection factor / varies between 0.1 and 1.0, meaning that 
10% to 100% of the available magnetic stress is released by reconnection. The total number 
of strands, Nm, goes from 20 to 1000, which is roughly the number of elemental magnetic 
flux tubes believed to be present in a coronal loop. We vary one parameter at a time keeping 
the rest of them at the intermediate values shown in bold face in Table [H 

As we stated above, the output of the CA model is an array of nanofiare energies 
dissipated in each strand as a function of time, with time steps of ~ 10^ s. Figure H] helps 
to illustrate the typical evolution of the system. For the example we use the parameter 
combination shown in bold face in Tableland a total duration of 72 timesteps (~ 20 hr). 
The upper panel shows the evolution of the mean energy of the system, defined as the sum 
of the magnetic energies of all the strands {Bf /Stc) divided by the number of strands. The 
middle panel shows the evolution of the mean nanofiare energy, which is the sum of the 
energies of all nanoflares occurred at each timestep divided by the number of strands. The 
evolution is such that the number of reconnection events increases from zero (rise phase) 




(7) 



Note that the parameter / affects the time required for a pair of adjacent strands to 
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to a quasi-stationary level (main phase). In the example of Figure H] this transition occurs 
around time step number 20. The system reaches this state abruptly, not asymptotically, and 
maintains it for as long as the driver follows the distribution of Figure |3]-b. Notice that once 
the system reaches the main phase the magnetic energy (Figure HI upper panel) maintains 
an approximately constant level, except for fluctuations due to the nanofiares. As soon as 
the driver begins to evolve towards the distribution of Figure [3]-c, the number of events 
decays until they almost disappear at the end of the system's evolution (Figure HJ middle 
panel). Likewise, once the decay phase begins the fluctuations of the system energy begin 
to decrease rapidly and stop. Summarizing, the evolution of the full system qualitatively 
follows the evolution of the observed coronal loops in terms of rise, main, and decay phases. 



3.3. Plasma response and light curve construction 

Since our main motivation is to compare the output of the CA model with the loops 
studied in Lopez Fuentes et al. (2007), the next step is to transform this output into an 
observed signal such as the light curve of a coronal loop. To do this, we simulate the plasma 
response using the EBTEL (Enthalpy-Based Thermal Evolution of Loops) model, developed 
by Klimchuk, Patsourakos & Cargill (2008). The temperature and density of the coronal 
plasma tend to be reasonably uniform along the magnetic field, and EBTEL computes the 
evolution of the spatially-averaged values. This approximate OD solution is quite similar to 
the exact solution given by ID hydrodynamic codes, but it takes 3-4 orders of magnitude 
less time to compute. EBTEL is therefore ideally suited for our study in which we simulate 
the evolution of thousands of individual strands. 

The main input of EBTEL is the heating rate profile, i.e., the variation of the heating 
rate with time in the strand. Here, we assume that nanofiares have a triangular profile of 
a given duration r. The upper panel of Figure [5] shows a series of "triangular" nanofiares 
that heat one of the strands of the CA model during a portion of the main phase. The time 
integral of each triangular heat function corresponds to the total energy provided by the 
particular nanoflare. The duration of the nanofiares shown in Figure [5] is 1000 s. This means 
that, in this case, the available heat is distributed (with a triangular profile) along a full time 
step of the CA model. The other two panels of Figure [5] show the output of EBTEL: the 
plasma temperature and density. During the time interval covered in the panels there are 
both isolated events (like the nanoflare starting at t = 4.3 x 10^ s, gray arrow), and "trains" 
of events occurring in rapid succession (such as the one starting at t = 4.7 x 10^ s, black 
arrow). Notice that while the plasma relaxes almost fully after the isolated event, during the 
succession of multiple events the plasma is still in a state of high density and temperature 
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due to the previous nanoflare when the next event begins. As we will see in the following 
sections, this has consequences in the way the plasma emits and the observed loop intensity 
evolves. 

The loops studied in Lopez Fuentes et al. (2007) have measured densities in the range 
0.9-1.9x10^ cm~'^, and measured temperatures between 1.2 and 2.1 MK. These values are 
reproduced by parameter values in the ranges given in Table [H 

Once we have the temperature and density evolution of each strand, we compute the 
observed intensity using the emission measure and the known temperature dependent sen- 
sitivity of the instrument. The emission measure per pixel of a single unresolved strand is 
EM = n'^AV, where AV is the volume of the strand portion covered by a pixel. This is 
roughly AV = Ipi^Ag, where Ipix is the pixel dimension and Ag is the cross sectional area of 
the strand. The contribution to the intensity in a pixel is then I = EM S{T), where S{T) 
is the SXI response function. We sum the contributions of all the strands to determine the 
intensity of the loop and, finally, we multiply this intensity by a factor that accounts for the 
proportion of the loop covered by a single pixel. The final output is the observed intensity 
of the loop per pixel per second as a function of time. The data are sampled using the 
instrument temporal resolution to simulate observed light curves. In Figure HI lower panel, 
we show a light curve obtained in this way for the model example of the upper and middle 
panels. Comparing the lower and middle panels it can be noticed that the light curve follows 
(though smoothed) the general fluctuations produced by the summed contribution of the 
nanoflares. 



4. Results 

The lower panels of Figure [1] show two model light curves obtained in the way described 
in Section 13. 3[ The cases shown have the same properties as the observed light curves of 
the upper panels. Given the random nature of the nanoflares, we do not expect an exact fit 
of the model to the observed data. For the comparison we instead consider the light-curve 
properties studied in Lopez Fuentes et al. (2007) (see Table [2]). 

The values of the model parameters used to construct the synthetic light curves of 
Figured] are, for Loop 1: = 37 G, tan^^ = 0.15 and / = 0.2; and for Loop 4: By = 20 G, 
tan 6'c = 0.2 and / = 0.1. For both cases we use = 100, and r = 1000 s. For the loop 
length and diameter we use the measures obtained in Lopez Fuentes et al. (2007): L = 110 
Mm for loop 1, L = 70 Mm for loop 4, and a diameter of 1.8 x 10^ cm for both. 

We compute for both observed and synthetic light curves the mean intensity of the main 
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phase (Im) and the duration and timescales of the rise and decay phases. To obtain these 
values we first identify, by visual inspection, the approximate start and end times of each 
phase. Using the same definitions as in Lopez Fuentes et al. (2007), we have for the rise 
phase timescale: 



Here, / is the mean intensity during the rise phase, AJ is the intensity variation and Atrise is 
the duration. A similar definition is used for the decay phase. As a proxy for the amplitude 
of the intensity fluctuations, we compute the root-mean-square (RMS) of the intensity 
variation during the main phase relative to its ten-point running average. In Table [2] we 
compare these properties for observed and modeled light curves. It can be seen that the 
model reproduces the light curve properties consistently. 

We remind the reader that the observed light curves shown in Figured] have been further 
processed (residual background removed, see Section [2]) from the versions studied in Lopez 
Fuentes et al. (2007). Therefore, the properties listed in Table [2] do not necessarily coincide 
with the corresponding values presented in the previous paper. 

Regarding the durations and timescales it is interesting to see that, since the light 
curves have initial intensities very close to zero, I in the Equation [8] is approximately A J/2. 
Therefore, 



Trise ~ „ • (9) 



We then expect the timescales to be approximately one half of the durations (compare the 
corresponding values in Table [2] for both the rise and decay phases). The relation from 
Equation is not valid, in general, for the light curves studied in Lopez Fuentes et al. 
(2007), in which there is still a remnant background component of the intensity. 

To investigate how the model parameters determine the properties of the loops, we vary 
each of them separately according to Table [H while keeping the rest fixed at the values shown 
in bold face (see Section [3]). In Figure Owe show log-log plots of the the mean energy per 
unit volume released in reconnection events, AE (see Equation [7]), versus different model 
parameters. The lines are least-square fits of the data, and their slopes are indicated in the 
panels. According to these plots, in our CA model the energy of the reconnection events 
approximately scales with -8^'°^ and (tan^c)^'^^- This can be easily explained by noting that 
in Equation [TJ when reconnection occurs, ABh ^ Be = B^tanOc- Therefore, it is expected 
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that AE oc 52 (tan6'c)2. 

The horizontal field decreases by an amount fABh/2 during a reconnection event. Sev- 
eral timesteps are typically required to rebuild the field to the critical level whereupon 
another event occurs. From equation (|2]) we see that the number of recovery timesteps is 
given by 



From this we obtain the mean nanofiare heating rate (averaged both in time and along the 
strand) : 



where tg is the time step duration and v = d/tg is the photospheric driving velocity. We have 
assumed that there is one reconnection event per nanofiare, whereas often there are two, but 
this is not important for identifying scaling relationships. In Figure [7] we show log-log plots 
of the mean heating rate versus By, tan 6'c, and L. The dependence on the parameters is: 
L~^-^, and (tan6'c)°'^^, which is clearly explained by Equation [TTl It is easy to see that 
(Q) does not depend on the number of strands, N^, since both the driver and the critical 
condition are applied to the strands individually. 

We can also analyze scaling relations for the properties of the light curves obtained 
with the model. For each combination of the parameters from Table [T] we obtain synthetic 
light curves following the procedure described in Section 13.31 In Figure M we present log-log 
plots of the mean intensity of the main phase, Im, as a function of the same parameters 
as in Figure [71 The scalings found in this case are: -B^'^^, and (tan^c)^'^^- Clearly, 

Im is related to the nanofiare energy, but the dependence is complex and is affected by 
both the detailed response of the plasma and the nonuniform temperature sensitivity of the 
instrument. is not linearly proportional to the total radiation emitted by the strand. 
To understand the origin of the dependence, let us now consider the case of quasi-static 
equilibrium. This is a reasonable representation if the time interval between nanoflares is 
short compared to a cooling time. When such conditions apply, the conductive and the 
radiation terms of the energy balance equation are comparable in magnitude to each other 
and to the mean heating rate, (Q) (Vesecky et al. 1979): 



A^r ~ — 7 tan 6'c. 



(10) 




(11) 



(Q) ^ -Ko 



2 



(12) 
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Here, is the coefficient of thermal conductivity, n is the electron number density, and Aq 
and h are constants. From the first part of Equation [12] we can obtain for the temperature: 

T oc {Qf L^" oc BfL^'\i^n e,f/\ (13) 

For the temperature range of the SXI loops studied in Lopez Fuentes et al. (2007), a rea- 
sonable single power value for the radiation loss is 6 ~ —0.5 (RTV, 1978). Using this, and 
equations [12] and [13] we have: 

^ oc ^ oc B,8/^L-3/7(tan 9,)^". (14) 

The observed loop intensity is / = n?S{T), where S{T) is the response function of 
the instrument, which we can express as S(T) = SoT"-. We found that a ~ 0.7 for the 
temperature range of the loops studied in Lopez Fuentes et al. (2007). Replacing equations [13] 
and [H] in the expression for the loop intensity, we obtain finally: 



/^oc5f9L-0-66(tane,)i-3^ (15) 

This expression is consistent with the scalings in Figure [8] Differences are to be expected 
since the strands are sometimes far from quasi-static equilibrium (Fig. [5]). 

It might be expected that the intensity of the main phase depends on A^^, but this is not 
the case. The reason is as follows. The intensity contribution of each strand is proportional 
to the volume of the strand subtended by the instrument pixel, which is proportional to the 
strand cross-section. Since we compute the strand cross section as the observed cross-section 
of the loop (1.8 X 10^ cm) divided by A^^, an increase of A^^ nieans a proportional decrease 
of each strand intensity, so the summed intensity of all strands remains unchanged. 

As we mentioned above, another important property of the light curves is the amplitude 
of the intensity fluctuations, which we quantify from the RMS of the intensity variation 
during the main phase. These fluctuations are related to the frequency and relative sizes of 
the nanofiares. As we briefly discussed in Section [n]the factor /, which relates to the fraction 
of the strand length that is involved in the reconnection, determines the waiting time for 
two neighbor strands to recover the critical condition. If / is small the waiting time will be 
short and the nanofiares will tend to be smaller and more frequent. On the other hand, as / 
tends to 1 the nanofiares will be larger and less frequent. Therefore, it is expected that the 
intensity fluctuations increase with /. The upper panel of Figure [1] shows precisely that. 
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The intensity fluctuations are also affected by the number of strands, Nm- Since the 
intensity of the loop is the summed contribution of all the strands, given the random nature 
of the release events, the larger the number of strands, the smoother the signal. This is 
shown in Figure IHl middle panel. We see that the RMS intensity fluctuation decreases as 
the number of strands increases. The effect levels off once a sufficiently large number of 
strands are present. 

We also explore how the nanofiare duration, r, affects the intensity fluctuations. The 
different values used are presented in Table [H last row. The results are shown in Figure [HI 
lower panel, where it can be seen that the RMS fluctuation also decreases as r increases. 
The reason is that individual strands evolve more slowly as the nanoflares get longer. For 
long enough nanoflares the strands behave in a quasi-static manner. As is the case with Nm, 
the effect levels off once the nanoflares are sufficiently long. 

Finally, although the dependence is much weaker than for the previous parameters, the 
RMS variation tends to increase with all the parameter changes that increase the intensity, 
i.e., the increase of and tan6'c, and the decrease of L. This can be understood in terms of 
the time it takes the plasma to cool fully after a nanofiare. The cooling time varies inversely 
with the nanofiare energy and directly with the loop length (Cargill & Klimchuk 2004). 
Shorter cooling times result in more bursty emission and greater fluctuation amplitudes. 
The increase of the fluctuation amplitude with the increasing intensity is consistent with 
the observational results obtained recently by Sakamoto et al. (2008, and references therein) 
and discussed also in Vekstein (2009). They also found that the amphtude of the intensity 
fluctuations of Yohkoh/SXT loops is approximately 8% of the value of the loop intensity. 
We find similar values in our own observations and the results obtained with the CA model 
(see the main phase RMS variation in Table |2]) . 

Another important property of the light curves is the rise duration, Atrise, which is the 
time it takes to reach the main phase after the intensity first starts to rise. Note that the 
intensity does not increase immediately when footpoint driving begins (see Figure HI middle 
and lower panels), because it takes time for the strands to reach the critical condition. The 
first strands to reach it are those that by chance have larger-than-average step sizes. We 
define the waiting time to be the elapsed time between the start of the evolution (t = 0) 
and the main phase beginning. It can be easily seen that is simply the number of steps 
required to reach the critical condition, \Bh\ = Bc/2, where the factor of 2 is because odd 
and even strands are tilted in opposite directions. The average \SBh\ per timestep that is 
injected by the driving is only half the value given by equation 1^, because one in four steps 
cause the stresses to decrease (Fig. [3]). Noting that Be = B^tsnaOc, we find that the waiting 
time is given by 
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in units of the timestep duration (1000 s). In Figure [TU] we show log-log plots of versus 
the parameters. The scalings are: L^-^^ and (tan6'c)°''^^. This is generally consistent with 
equation f|T6|) . and we suggest that small deviations are due to the subjective nature of 
identifying the start of the main phase. Notice that is expected to be independent of 
By, which we have confirmed to be true from our simulations. We emphasize again that the 
waiting time is different from the rise duration. The rise duration is more closely related to 
the spread of the driver distribution (e.g., the Gaussian width in Fig. 13b) than it is to the 
time required for the average strand to reach critical. 

The duration of the decay phase depends on both the form of the final driver distribution 
(Fig. Eb) and the timescale for transitioning from the initial to final distributions. We leave 
the details for a later study. 



5. Discussion 

We develop a pseudo-ID model based on a cellular automaton (CA) that reproduces 
the main properties of the evolution of coronal loops observed in soft X-rays. As we describe 
in Section El the model is based on the idea that coronal heating is due to the sudden recon- 
nection of tangled magnetic strands (Parker 1988). Using typical solar ranges for the model 
parameters we qualitatively reproduce the observed intensity amplitude, fluctuations, and 
evolution timescales of soft X-ray loops studied in Lopez Fuentes et al. (2007). The syn- 
thetic loops obtained with the CA-I-EBTEL model combination also reproduce the measured 
physical properties of the observed loop plasma such as density and temperature. 

Although our model is related to usual models of self-organized criticality (SOC, see e.g. 
Charbonneau et al. 2001), it differs in some important aspects. In most SOC models the 
rules for the critical condition and the redistribution of the field among neighboring strands 
are based on mean values of the field strength. The strength of each strand is compared 
with the mean value of its neighbors, and if the critical condition is fulfilled the field is 
redistributed equally among those neighbors. One may classify this type of model isotropic. 
As we describe in Section [31 our CA model is anisotropic in this respect. The other important 
difference is that SOC models use null boundary conditions. The mesh points (strands) at 
the boundary are set to have zero magnetic field (what we call Bh in Section ^ at all times. 
Here, we use instead periodic boundary conditions. It may be partly due to these differences 
that the nanofiares produced with our model do not follow power law distributions, which 
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are characteristic of SOC models (for an anisotropic SOC model associated with power law 
distributions, see e.g. Morales & Charbonneau 2008). We must stress that our main goal 
here is not to obtain such distributions, but instead to explain the evolution of loops in the 
frame of the nanoflare model. It is worth pointing out, though, that our CA model reaches 
a steady critical state in much the same way as SOC models do (see Figure H]). 

As we mention in Section [1], in recent years there has been a debate over the issue 
of whether coronal loops are isothermal or multithermal. The argument goes that if loops 
are isothermal they are likely to be monolithic structures, while multithermality implies that 
loops are composed of different unresolved strands evolving independently. There seems to be 
evidence in favor of both kinds of observations (see e.g., Schmelz et al. 2009), so a consensus 
is forming that both scenarios actually occur. Let us analyze this possibility in terms of the 
model presented here. For instance, if the reconnection mechanism is not very efficient (e.g., 
a small / factor, see Section [3]) and all nanoflares are approximately the same size, then 
the heating injection will be difficult to distinguish from a single source of constant heating. 
In that case, the strand will keep more or less the same temperature along the main phase 
of the evolution. If all the strands that comprise the loop have similar physical conditions 
(it is expected that they all have approximately the same field strength, for example), then 
the loop will be nearly isothermal. These loops can be easily taken as examples of quasi- 
static evolution. On the other hand, if nanoflares occur less frequently, so the waiting time 
between consecutive nanoflares is not negligible, the strands have time to cool to lower 
temperatures between one nanoflare and the next. Since nanoflares occur randomly, at any 
given time during the evolution of the loop there will be strands simultaneously having 
different temperatures. That would correspond to a multithermal loop observation. 

There is yet a third possibility. If the nanoflares occur only once in each strand, with- 
out repeating, and if they all occur at approximately the same time, then the loop will 
appear roughly isothermal at any given moment. The temperature will decrease on a cool- 
ing timescale, however, and the loop will be short lived. We call this scenario a short duration 
nanoflare "storm" (Klimchuk 2009). Many warm EUV loops can be explained in this way. 
Hot soft X-ray loops tend to live much longer than a cooling time, on the other hand (Lopez 
Fuentes et al. 2007). If they are indeed isothermal, then they can only be explained by 
rapidly repeating nanoflares or by truly steady heating. 

The three alternative scenarios described above relate to the question weather the same 
loop is observable by instruments that are sensitive to different temperatures. In the case of 
nanoflares that repeat rapidly in each strand of the loop bundle, the loop will be visible only 
over a narrow range of temperatures. It will be seen either in hot soft X-ray emission or in 
warm EUV emission, but not both. In the case of nanoflares that repeat with a time delay 



- 18 - 



longer than the coohng time, the loop will be simultaneously visible in both emissions. In 
the final case of a short duration nanoflare storm, the loop will be visible first in soft X-rays 
and then in the EUV. 

Observational evidence has been presented for all three possibilities. The question of 
which is most prevalent is not yet answered. Loops produced by short duration storms 
have been clearly identified and can be found at many locations within active regions (e.g., 
Ugarte-Urra et al. 2006, 2009). On the other hand, soft X-ray loops are most pronounced in 
the central cores of active regions, where EUV loops are less common (Schmieder et al. 2004). 
There are also reports of soft X-ray and EUV loops that are almost but not exactly co-spatial 
(Nitta 2000, Nagata et al. 2003). We attempted to study the warm emission associated with 
the SXI loops studied in Lopez Fuentes et al. (2007). Unfortunately, TRACE data are not 
available, so we had to rely on lower cadence observations from SOHO/EIT. Since only 
four images per day were taken in each channel, we could only compare the hot and warm 
emission at a few times during the loop evolution. We found that some of the SXI loops 
had an EIT counterpart but others did not. Clearly, more systematic observations of loops 
in different wavelengths are needed. 

It is possible that there are two well differentiated populations of AR loops. Hotter 
and longer-lived loops are preferentially located at the core of ARs (Antiochos et al. 2003, 
Ugarte-Urra et al. 2009), while the periphery is dominated by longer, shorter- lived, cooler 
loops (Del Zanna & Mason 2003). If this is indeed the case in general, then the different 
locations and evolutions of these loop classes may be related to different physical properties 
at their origin. The model we present in this paper applies specifically to hot loops. Whether 
it can explain warm loops at the periphery of ARs is unclear. We can nonetheless speculate 
on the difference in these two regions. The strength of the coronal magnetic field is generally 
greater in the core of an AR and this will tend to produce stronger heating. Our model 
also depends fundamentally on the properties of footpoint shuffling. It may be that these 
properties vary systematically across the active region. This is an important question that 
we hope will be answered with the latest magnetogram observations from SOT/Hinode. 
These observations will hopefully also confirm our proposal that the slow decay phase of hot 
loops is associated with the dispersal of magnetic footpoints to the point where magnetic 
field tangling is no longer efficient. 

In Section H] we found a series of scaling laws relating observed properties of the loop 
evolution to different physical parameters of the model. Some of these parameters, such 
as the loop length, can actually be obtained directly from the observations. The strength 
and tilt of the magnetic field at the coronal base can be inferred by combining photospheric 
magnetogram observations with models of the coronal structure (see e.g., Lopez Fuentes et 
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al. 2006, Mandrini et al. 2000). Therefore, these scahng laws are predictions of the model 
that can be tested observationally. 

Prom the discussion in the previous paragraphs we conclude that among future impor- 
tant investigations is the systematic study of loop observations in different wavelengths and 
the analysis of possible scalings of the evolution parameters (timescales, mean intensities, 
fluctuations) with the loop physical properties such as length and magnetic field. From the 
theoretical side, the model presented here is still very simple, and more sophisticated ver- 
sions should be developed, such as 2D models that better account for the vector nature of 
the magnetic field and the fact that each magnetic strand is wrapped around multiple other 
strands. We have already started work in this direction. 

The authors wish to thank Paul Charbonncau for enriching discussions and the anony- 
mous referee for useful comments and suggestions. This work was supported by the NASA 
Supporting Research and Technology and Living With a Star Programs. 
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Table 1: Numerical values used for the parameters of the model. 
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Table 2: Compared properties of the observed and modeled lightcurves shown if Figured] 



Loop 1 Loop 4 

Obs. Model Obs^ Model 



Rise phase duration (hr) 


3.35 


3.65 


3.74 


3.16 


Rise phase timescale (hr) 


1.97 


1.88 


2.17 


1.75 


Main phase intensity (DN/pix.sec) 


6.82 


7.59 


2.99 


2.82 


Main phase RMS 


0.10 


0.10 


0.13 


0.09 


Decay phase duration (hr) 


2.98 


3.99 


4.51 


4.19 


Decay phase timescale (hr) 


-1.51 


-2.20 


-2.47 


-2.4J 
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Fig. 1. — Evolution of observed and synthetic X-ray loops. The light curves obtained with 
the CA model (lower panels) reproduce the main properties of the evolution of the observed 
loops (upper panels). See Section HI 
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Fig. 2. — Schematic description of the CA model, (a) Initially, all strands are parallel and 
there is no magnetic stress between neighbors. The horizontal planes represent two distant 
portions of the photosphere connected by the strands, (b) Photospheric motions displace the 
strand footpoints increasing the mutual inclinations and producing magnetic stress, (c) We 
associate the inclination of the strands (6) with the appearence of a horizontal component 
of the magnetic field (see Section [XT]) . 
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Fig. 3. — Panel (a): the gray arrow indicates the footpoint displacement from time step (t—l) 
to time step t and the black arrows indicate the possible further evolution from timestep 
t to timestep (t + 1), which may increase (3/4 probability) or cancel (1/4 probability) the 
displacement of the previous step. The driver distribution of panel (b) is consistent with 
these probabilities. To simulate the footpoint dispersion, the distribution of panel (b) slowly 
evolves to the distribution of panel (c) (see Section 13. 2p 
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Fig. 4. — Example of the CA model evolution. Upper panel: mean magnetic energy (per 
strand) versus time step number. Middle panel: idem for the mean nanoflare energy (per 
strand, see Section 13. 2p . Lower panel: light curve obtained from the model of the upper 
panels (see Section 13. 3p . The 3 panels have a precise vertical alignment, so that 20 hr 
corresponds to time step 72. 
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Fig. 5. — Plasma response to a series of nanoflares modeled with the EBTEL code. The gray 
arrow indicates the start time of a isolated nanoflare, while the black arrow corresponds to 
the start of a "train" of events occurring in a rapid succession (see Section 1331) . 
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Fig. 6. — Log-log plots of the mean reconnection event energy, AE, versus relevant param- 
eters of the model (see Section Hj). 
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Fig. 7. — Log-log plots of the mean nanoflare heating rate, {Q), versus parameters of the 
model (see Section H]). 
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Fig. 8. — Log-log plots of the main phase mean intensity, Im, versus relevant parameters of 
the model (see Section H]). 
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Fig. 9. — Scatter plots of the root mean square variation of the main phase intensity, RMS, 
versus relevant parameters of the model. 
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Fig. 10. — Log-log plots of the waiting time, (see Section Hj), versus relevant parameters 
of the model. 



